BIC metabolomics data analysis

In this document we present the joint analysis of the PASS1A metabolomics datasets.

Load all datasets and setup the session

Load the data from the cloud, including: phenotypic data, metabolomic datasets, and metabolomics dictionary.

source("~/Desktop/repos/motrpac-bic-norm-qc/tools/supervised_normalization_functions.R")
source("~/Desktop/repos/motrpac-bic-norm-qc/tools/unsupervised_normalization_functions.R")
source("~/Desktop/repos/motrpac-bic-norm-qc/tools/gcp_functions.R")
source("~/Desktop/repos/motrpac-bic-norm-qc/tools/association_analysis_methods.R")
source("~/Desktop/repos/motrpac-bic-norm-qc/tools/data_aux_functions.R")
source("~/Desktop/repos/motrpac/tools/prediction_ml_tools.R")
library(randomForest) # for classification tests
# Load the dmaqc data
merged_dmaqc_data =  load_from_bucket("merged_dmaqc_data2019-10-15.RData",
    "gs://bic_data_analysis/pass1a/pheno_dmaqc/",F)

merged_dmaqc_data = merged_dmaqc_data[[1]]
rownames(merged_dmaqc_data) = as.character(merged_dmaqc_data$vial_label)
# define the tissue variable
merged_dmaqc_data$tissue = merged_dmaqc_data$sampletypedescription
# define the time to freeze variable
merged_dmaqc_data$time_to_freeze = merged_dmaqc_data$calculated.variables.time_death_to_collect_min + 
  merged_dmaqc_data$calculated.variables.time_collect_to_freeze_min
# Add blood freeze times
blood_samples = 
  merged_dmaqc_data$specimen.processing.sampletypedescription ==
  "EDTA Plasma"
blood_freeze_time = 
  as.difftime(merged_dmaqc_data$specimen.processing.t_freeze,units = "mins") -
  as.difftime(merged_dmaqc_data$specimen.processing.t_edtaspin,units="mins")
blood_freeze_time = as.numeric(blood_freeze_time)
merged_dmaqc_data$time_to_freeze[blood_samples] = 
  blood_freeze_time[blood_samples]
# Load our parsed metabolomics datasets
metabolomics_bucket_obj = load_from_bucket(
  file = "metabolomics_parsed_datasets_pass1a_external1.RData",
  bucket = "gs://bic_data_analysis/pass1a/metabolomics/")

metabolomics_parsed_datasets = metabolomics_bucket_obj[[1]]
# Load the processed, merged, normalized datasets
metabolomics_processed_datasets = load_from_bucket(
  file="metabolomics_processed_datasets11182019.RData",
  bucket = "gs://bic_data_analysis/pass1a/metabolomics/"
)[[1]]

# Define the variables to be adjusted for:
biospec_cols = c(
  "acute.test.distance",
  "calculated.variables.time_to_freeze",
  # "calculated.variables.edta_coll_time", # no need - see code above for blood
  "bid" # required for matching datasets
  )
differential_analysis_cols = c(
  "animal.registration.sex",
  "animal.key.timepoint",
  "animal.key.is_control"
)
pipeline_qc_cols = c("sample_order")

Log-transform: effect on variance

Some sites do not use the log transformation on their dataset. In this section we plot the coefficient of variation as a function of the mean instensity. We take a single dataset as an example to show how log-transformed data have reduced dependency and smoother plots.

As an additional analysis we also plot the number of missing values per metabolite as a function of its mean intensity. We show that while there is high correlation some missing values appear in fairely high intensities. This is important for imputation as some sites use some fixed low value instead of knn imputation.

# Plot cv vs means
library(gplots)
d = metabolomics_parsed_datasets[["white_adipose_powder,metab_u_hilicpos,unnamed"]]
dx = d$sample_data
CoV<-function(x){return(sd(x,na.rm = T)/mean(x,na.rm=T))}
dmeans = apply(dx,1,mean,na.rm=T)
CoVs = apply(dx,1,CoV)
inds = !is.na(CoVs)
df = data.frame(Mean_intensity = dmeans[inds],CoV = CoVs[inds])
plot(CoV~Mean_intensity,df,cex=0.5,pch=20,main="Raw data")
lines(lowess(CoV~Mean_intensity,df),lty=2,lwd=2,col="blue")

# Repeat after log2
dx = log(1+d$sample_data,base=2)
dmeans = apply(dx,1,mean,na.rm=T)
CoVs = apply(dx,1,CoV)
inds = !is.na(CoVs)
df = data.frame(Mean_intensity = dmeans[inds],CoV = CoVs[inds])
plot(CoV~Mean_intensity,df,cex=0.5,pch=20,main="Log2 data")
lines(lowess(CoV~Mean_intensity,df),lty=2,lwd=2,col="blue")

# Plot number of NAs vs intensity mean
dx = log(1+d$sample_data,base=2)
dmeans = apply(dx,1,mean,na.rm=T)
num_nas = rowSums(is.na(dx))
df = data.frame(Num_NAs = num_nas[inds],Mean_intensity = dmeans[inds])
rho = cor(df$Num_NAs,df$Mean_intensity,method="spearman")
rho = format(rho,digits=2)
plot(Num_NAs~Mean_intensity,df,cex=0.5,pch=20,
     main=paste("Spearman:",rho))

Log-transform: effect on differential analysis

Correct the meta-data, for the differential analysis:

for(datasetname in names(metabolomics_processed_datasets)){
  curr_data = metabolomics_processed_datasets[[datasetname]]
  # organize the metadata
  curr_meta = merged_dmaqc_data[colnames(curr_data),
        union(biospec_cols,differential_analysis_cols)]
  # remove metadata variables with too many NAs
  na_counts = apply(is.na(curr_meta),2,sum)
  curr_meta = curr_meta[,na_counts/nrow(curr_meta) < 0.1]
}

Error in `[.data.frame`(curr_meta, , na_counts/nrow(curr_meta) < 0.1) : 
  undefined columns selected

Untargeted data are typically log-transformed and analyzed using linear models. On the other hand, concentration data are sometimes analyzed with the same type of models but using the original data. This raises a problem if we wish to compare exact statistics from these data. In this section we perform residual analysis for single metabolites. Our goal is to identify if concentration data behaves “normally” when not log-transformed. The analysis below examines the residuals of the data after fitting linear models for each metabolite, adjusting for freeze time and sex. We then compare the results with and without the log-transformation, counting the number of metabolites with a significant evidence for non-normally distributed residuals.

Trageted vs. Untargeted: single metabolite comparison

Compute statistics for each dataset

Compare overlaps, effect sizes, and correlations within tissues. Compare targeted-untargeted pairs only. For differential analysis we use the same model as in the analysis above.

We next transform the data above into tables that contain data for each combination of metabolite, time point, and tissue. These are then used for different meta-analyses: (1) a simple random effects analysis, (2) random effects with a binary covariate indicating if a dataset is targeted or untargeted, (3) redo the RE model of (1) with the targeted data only, and (4) redo the RE model of (1) with the untargeted data only.

for(tissue in unique(single_metabolite_de$tissue)){
  for(tp in unique(single_metabolite_de$tp)){
    curr_subset = single_metabolite_de[
      single_metabolite_de$tissue==tissue &
        single_metabolite_de$tp==tp,]
    for(metabolite in unique(curr_subset$metabolite)){
      curr_met_data = curr_subset[
        curr_subset$metabolite==metabolite,]
      curr_met_data$var = curr_met_data$Std^2
      re_model1 = NULL;re_model2=NULL
      re_model_tar = NULL;re_model_untar = NULL
      try({re_model1 = rma.uni(curr_met_data$Est,curr_met_data$var)})
      try({re_model2 = rma.mv(curr_met_data$Est,curr_met_data$var,
                         mods=curr_met_data$is_targeted)})
      try({re_model_tar = rma.uni(
        curr_met_data[curr_met_data$is_targeted==1,"Est"],
        curr_met_data[curr_met_data$is_targeted==1,"var"]
      )})
      try({re_model_untar = rma.uni(
        curr_met_data[curr_met_data$is_targeted==0,"Est"],
        curr_met_data[curr_met_data$is_targeted==0,"var"]
      )})
      meta_analysis_stats[[paste(metabolite,tissue,tp,sep=",")]] = 
        list(curr_met_data=curr_met_data,re_model1=re_model1,
            re_model2 = re_model2,re_model_tar=re_model_tar,
            re_model_untar = re_model_untar)
    }
  }
}

Error in rma.uni(curr_met_data$Est, curr_met_data$var) : 
  Fisher scoring algorithm did not converge. See 'help(rma)' for possible remedies.

Comparison results

We now show some plots to summarize the comparison.

Dataset coverage

We first plot the number and percentage of metabolites in the targeted datasets that are measured in at least one untargeted dataset.

Spearman correlations

We examine the average absolute correlation between the platforms (within tissues). Whenever two platforms share more than a single metabolite we plot both the average correlation between the same metabolites and between other metabolites. Adding the average correlation between platforms but with different metabolites is important as it gives some perspective to what a significant correlation is. That is, in many cases below, the average correlation may be greater than expected.

Meta-analysis of differential effects

Here we go into the single metabolites comparison in greater detail (again, within tissues).

We start with a few examples. Here are the results for lactate in plasma.

lact_res = meta_analysis_stats[
  grepl("lact",names(meta_analysis_stats),ignore.case = T) &
    grepl("plasma",names(meta_analysis_stats),ignore.case = T)
]
lact_res_hours = sapply(names(lact_res),
            function(x)as.numeric(strsplit(x,split=",")[[1]][3]))
lact_res = lact_res[order(lact_res_hours)]
for(lact_example in names(lact_res)[1:6]){
  curr_labels = gsub("plasma,","",
                     meta_analysis_stats[[lact_example]][[1]][,1])
  forest(meta_analysis_stats[[lact_example]]$re_model1,
       slab = curr_labels,
       main = lact_example,xlab = "Log fc",
       col = "blue",cex = 1.1)
}

We can now check the same analysis for liver:

lact_res = meta_analysis_stats[
  grepl("lact",names(meta_analysis_stats),ignore.case = T) &
    grepl("liver",names(meta_analysis_stats),ignore.case = T)
]
lact_res_hours = sapply(names(lact_res),
            function(x)as.numeric(strsplit(x,split=",")[[1]][3]))
lact_res = lact_res[order(lact_res_hours)]
for(lact_example in names(lact_res)[1:6]){
  curr_labels = gsub("liver_powder,","",
                     meta_analysis_stats[[lact_example]][[1]][,1])
  forest(meta_analysis_stats[[lact_example]]$re_model1,
       slab = curr_labels,
       main = lact_example,xlab = "Log fc",
       col = "blue",cex = 1.1)
}

We now move to a systematic comparison of the datasets. We start with a naive comparison, checking if metabolites are consistently being discovered using a simple \(p=0.001\) threshold for significance.

# Naive comparison
thr = 0.001
naive_analysis_tables = lapply(meta_analysis_stats,
    function(x)table(x$curr_met_data[,"Pvalue"]<thr,
                     x$curr_met_data[,"is_targeted"]))
table_with_sig_results = naive_analysis_tables[
  sapply(naive_analysis_tables,nrow)>1
]
percent_with_sig = 
  length(table_with_sig_results)/length(naive_analysis_tables)
sig_results = sapply(table_with_sig_results,
      function(x)paste(x["TRUE",c("0","1")]>0,collapse=","))
sig_results_counts = table(sig_results)
print("Counting the number of metabolites with p<0.001")

[1] "Counting the number of metabolites with p<0.001"

print(paste("Significant in targeted only:",sig_results_counts["FALSE,TRUE"]))

[1] "Significant in targeted only: 137"

print(paste("Significant in untargeted only:",sig_results_counts["TRUE,FALSE"]))

[1] "Significant in untargeted only: 190"

print(paste("Significant in both:",sig_results_counts["TRUE,TRUE"]))

[1] "Significant in both: 100"

As a more rigorous analysis, we use the meta-analyses to obtain useful statistics for comparison.

# Extract some useful statistics per analysis
# P-value for the difference between targeted and untargeted
targeted_diff_p = 
  sapply(meta_analysis_stats,function(x)x$re_model2$pval[2])
# P-values - targeted vs. untargeted
pvals_tar = sapply(meta_analysis_stats,function(x)x$re_model_tar$pval)
pvals_untar = sapply(meta_analysis_stats,function(x)x$re_model_untar$pval)
pvals_untar = unlist(pvals_untar[sapply(pvals_untar,length)>0])
significant_in = rep("None",length(pvals_untar))
significant_in[pvals_tar<0.001] = "Targeted"
significant_in[pvals_untar<0.001] = "Untargeted"
significant_in[pvals_tar<0.001 & pvals_untar<0.001] = "Both"
significant_diff = targeted_diff_p<0.001
df = data.frame(
  targeted = -log10(pvals_tar),
  untargeted = -log10(pvals_untar),
  significant_in = significant_in,
  significant_diff = significant_diff
)
rho = cor(pvals_tar,pvals_untar,method = "spearman")
rhop = cor.test(pvals_tar,pvals_untar,method = "spearman")$p.value
print(
  ggplot(df, aes(x=targeted, y=untargeted,
                 shape=significant_diff, color=significant_in)) +
    geom_point() +
    ggtitle(paste("-log10 p-values, spearman:",format(rho,digits=2),
                  "(p=",format(rhop,digits=3),")"))
)

print("### Summary of differences in RE models ###")

[1] "### Summary of differences in RE models ###"

print("Model is significant at p<0.001:")

[1] "Model is significant at p<0.001:"

print(table(df$significant_in))

      Both       None   Targeted Untargeted 
       169       2029        180        275 

print("Adding is_targeted as a covariate has p<0.001:")

[1] "Adding is_targeted as a covariate has p<0.001:"

print(table(df$significant_diff))

FALSE  TRUE 
 2362   291 

# Betas - targeted vs. untargeted
betas_tar = sapply(meta_analysis_stats,function(x)x$re_model_tar$beta[1,1])
betas_untar = sapply(meta_analysis_stats,function(x)x$re_model_untar$beta[1,1])
betas_untar = unlist(betas_untar[sapply(betas_untar,length)>0])
df = data.frame(
  targeted = betas_tar,
  untargeted = betas_untar,
  significant_in = significant_in,
  significant_diff = significant_diff
)
rho = cor(betas_untar,betas_tar,method = "spearman")
rhop = cor.test(betas_untar,betas_tar,method = "spearman")$p.value
print(
  ggplot(df, aes(x=targeted, y=untargeted,
                 shape=significant_diff, color=significant_in)) +
    geom_point() +
    ggtitle(paste("Effect sizes, spearman:",format(rho,digits=2),
                  "(p=",format(rhop,digits=3),")"))
)

From the plots above we take the most extreme examples and examine their forest plots.

Targeted vs. untargeted: comparison as a prediction task

Use 5-fold cross validation for analysis within tissues. For each pair of targeted and untargeted datasets from the same tissue, we use the untargeted data as the predictive features and all metabolites in the targeted datasets as the dependent variables. The code below uses feature selection and random forests to train the predictive models.

nfolds = 5
prediction_analysis_results = list()
for(nn1 in names(metabolomics_processed_datasets)){
  nn1_tissue = strsplit(nn1,split=",")[[1]][1]
  nn1_tissue = gsub("_powder","",nn1_tissue)
  if(grepl("untargeted",nn1)){next}
  for(nn2 in names(metabolomics_processed_datasets)){
    if(nn2 == nn1){next}
    if(!grepl("untargeted",nn2)){next}
    nn2_tissue = strsplit(nn2,split=",")[[1]][1]
    nn2_tissue = gsub("_powder","",nn2_tissue)
    nn2_dataset = strsplit(nn2,split=",")[[1]][2]
    if(nn1_tissue!=nn2_tissue){next}
    print(paste("features from:",nn2))
    print(paste("labels from:",nn1))
    # get the numeric datasets and their annotation
    y = metabolomics_processed_datasets[[nn1]]$sample_data
    x = metabolomics_processed_datasets[[nn2]]$sample_data
    # align the sample sets
    bid_y = merged_dmaqc_data[colnames(y),"bid"]
    bid_x = merged_dmaqc_data[colnames(x),"bid"]    
    # step 1: merge samples from the same BID
    if(length(unique(bid_x))!=length(bid_x)){
      x = aggregate_repeated_samples(x,bid_x)
    }
    else{
      colnames(x) = bid_x
    }
    if(length(unique(bid_y))!=length(bid_y)){
      y = aggregate_repeated_samples(y,bid_y)
    }else{
      colnames(y) = bid_y
    }
    # step 2: use the shared bio ids
    shared_bids = as.character(intersect(colnames(y),colnames(x)))
    x = t(as.matrix(x[,shared_bids]))
    y = t(as.matrix(y[,shared_bids]))
    # At this point x and y are over the same BIDs, now we add the metadata
    y_meta = unique(metabolomics_processed_datasets[[nn1]]$sample_meta_parsed)
    rownames(y_meta) = y_meta$bid
    y_meta = y_meta[shared_bids,]
    
    # take the covariates (ignore distances)
    curr_cov_cols = intersect(colnames(y_meta),biospec_cols[2])
    curr_covs = data.frame(y_meta[,curr_cov_cols])
    names(curr_covs) = curr_cov_cols
    curr_covs$sex = y_meta$animal.registration.sex # add sex
    # add the covariates into x
    x = cbind(x,curr_covs)
    
    # Run the regressions
    folds = sample(rep(1:nfolds,(1+nrow(x)/nfolds)))[1:nrow(x)]
    numFeatures = min(ncol(x),2000)
    preds = c();real=c()
    for(i in 1:ncol(y)){
      if( i %% 10 == 0){print(paste("analyzing metabolite number:",i))}
      y_i = y[,1]
      i_preds = c();i_real=c()
      for(j in 1:nfolds){
        tr_x = x[folds!=j,]
        tr_yi = y_i[folds!=j]
        te_x = x[folds==j,]
        te_y = y_i[folds==j]
        # random forest
        # model = randomForest(tr_yi,x=tr_x,ntree = 20)
        # te_preds = predict(model,newdata = te_x)
        model = feature_selection_wrapper(tr_x,tr_yi,
                   coeff_of_var,randomForest,
                   topK = numFeatures,ntree=50)
        te_preds = predict(model,newdata = te_x)
        i_preds = c(i_preds,te_preds)
        i_real = c(i_real,te_y)
      }
      preds = cbind(preds,i_preds)
      real = cbind(real,i_real)
    }
    colnames(preds) = colnames(y)
    colnames(reals) = colnames(y)
    currname = paste(nn1,nn2,sep=";")
    prediction_analysis_results[[currname]] = list(
      preds = preds,real=real
    )
  }
}
save_to_bucket(prediction_analysis_results,
               file="tar_vs_untar_prediction_analysis_results.RData",
               bucket = "gs://bic_data_analysis/pass1a/metabolomics/")

We now take the predicted and real values and estimate the prediction accuracy in different ways.

prediction_analysis_results =  load_from_bucket(
  "tar_vs_untar_prediction_analysis_results.RData",
    "gs://bic_data_analysis/pass1a/metabolomics/",F)
prediction_analysis_results = prediction_analysis_results[[1]]

results_metrics = list()
for(nn in names(prediction_analysis_results)){
  preds = prediction_analysis_results[[nn]]$preds
  real = prediction_analysis_results[[nn]]$real
  tar_name = strsplit(nn,split=";")[[1]][1]
  untar_name = strsplit(nn,split=";")[[1]][2]
  y = metabolomics_processed_datasets[[tar_name]]$sample_data
  colnames(preds) = rownames(y)
  colnames(real) = rownames(y)
  tar_name = simplify_metab_dataset_name(tar_name)
  untar_name = simplify_metab_dataset_name(untar_name)
  currtissue = strsplit(tar_name,split=",")[[1]][1]
  tar_name = gsub(paste(currtissue,",",sep=""),"",tar_name)
  untar_name = gsub(paste(currtissue,",",sep=""),"",untar_name)
  if(! currtissue %in% names(results_metrics)){
    results_metrics[[currtissue]] = list()
  }
  if(! tar_name %in% names(results_metrics[[currtissue]])){
    results_metrics[[currtissue]][[tar_name]] = list()
  }
  
  rhos = format(diag(cor(preds,real,method="spearman")),digits=3)
  rhos = as.numeric(rhos)
  SEs = colSums((preds-real)^2)
  MSEs = SEs / nrow(preds)
  RMSE = sqrt(MSEs)
  rMSE = MSEs / apply(y,1,var)
  CoVs = apply(y,1,sd) / apply(y,1,mean)
  discCoVs = cut(CoVs,breaks = 2,ordered_result = T)
  
  results_metrics[[currtissue]][[tar_name]][[untar_name]] = data.frame(
    rhos,MSEs,RMSE,rMSE,CoVs,discCoVs
  )
}

We now present a few summary plots.

for(tissue in names(results_metrics)){
  for(tar in names(results_metrics[[tissue]])){
    l = results_metrics[[tissue]][[tar]]
    rho_vs_cv = c()
    for(untar in names(l)){
      m = l[[untar]][,c("rhos","discCoVs")] # take the current matrix
      m = cbind(rep(untar,nrow(m)),m)
      m$discCoVs = as.numeric(m$discCoVs)
      rho_vs_cv = rbind(rho_vs_cv,m)
    }
    colnames(rho_vs_cv)[1] = "dataset"
    boxplot(rhos~discCoVs:dataset,data=rho_vs_cv,las=2,
            ylab="Spearman",xlab = "",ylim=c(0,1),
            main = paste(tissue,tar,sep=","))
  }
}

As additional references, we train below additional models. First, we check the prediction of naive models that use technical and clinical covariates only. Second, we use multi-task regression and deep learning models.

cov_prediction_analysis_results = list()
for(nn1 in names(metabolomics_processed_datasets)){
  nn1_tissue = strsplit(nn1,split=",")[[1]][1]
  nn1_tissue = gsub("_powder","",nn1_tissue)
  if(grepl("untargeted",nn1)){next}
  print(nn1)
  y = metabolomics_processed_datasets[[nn1]]$sample_data
  y_vials = colnames(y)
  bid_y = merged_dmaqc_data[colnames(y),"bid"]
  colnames(y) = bid_y
  y = t(as.matrix(y))
  if(ncol(y)>1000){next}
  cov_cols = c("animal.registration.sex",
             "acute.test.weight",
             "acute.test.distance",
             "animal.key.timepoint")
  covs = merged_dmaqc_data[y_vials,cov_cols]
  x = covs
  
  # Run the regressions
  folds = sample(rep(1:nfolds,(1+nrow(x)/nfolds)))[1:nrow(x)]
  numFeatures = min(ncol(x),2000)
  preds = c();real=c()
  for(i in 1:ncol(y)){
    y_i = y[,1]
    i_preds = c();i_real=c()
    for(j in 1:nfolds){
      print(j)
      tr_x = x[folds!=j,]
      tr_yi = y_i[folds!=j]
      te_x = x[folds==j,]
      te_y = y_i[folds==j]
      # random forest
      model = randomForest(tr_yi,x=tr_x,ntree = 20)
      te_preds = predict(model,newdata = te_x)
      i_preds = c(i_preds,te_preds)
      i_real = c(i_real,te_y)
    }
    preds = cbind(preds,i_preds)
    real = cbind(real,i_real)
  }
  cov_prediction_analysis_results[[nn1]] = list(
      preds = preds,real=real
    )
}

# preds = c();real=c()
# for(j in 1:nfolds){
#   tr_x = x[folds!=j,]
#   tr_y = y[folds!=j,]
#   te_x = x[folds==j,]
#   te_y = y[folds==j,]
#   model = MTL_wrapper(tr_x,tr_y,type="Regression", Regularization="L21")
#   te_preds = predict(model,te_x)
#   real = rbind(real,te_y)
#   preds = rbind(preds,te_preds)
# }
# diag(cor(preds,real))

# Using PLS regression
# library(pls)
# pls_model = plsr(y~x,ncomp = 5,validation="LOO")
# eval = MSEP(pls_model)
# 
# y_pca = prcomp(y)
# plot(y_pca)
# explained_var = y_pca$sdev^2/sum(y_pca$sdev^2)
# y_pca_matrix = y_pca$x[,1:10]
# 
# # regress out sex, weight
# 
# get_explained_variance_using_PCA(x,y)
# x = apply(x,2,regress_out,covs=covs)
# y = apply(y,2,regress_out,covs=covs)
# get_explained_variance_using_PCA(x,y)